Current-induced nonequilibrium vibrations in single-molecule devices 
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Finite-bias electron transport through single molecules generally induces nonequilibrium molecular 
vibrations (phonons). By a mapping to a Fokker-Planck equation, we obtain analytical scaling forms 
for the nonequilibrium phonon distribution in the limit of weak electron-phonon coupling A within 
a minimal model. Remarkably, the width of the phonon distribution diverges as ~ X~ a when the 
coupling decreases, with voltage-dependent, non-integer exponents a. This implies a breakdown 
of perturbation theory in the electron-phonon coupling for fully developed nonequilibrium. We 
also discuss possible experimental implications of this result such as current-induced dissociation of 
molecules. 
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I. INTRODUCTION 

The vision of molecular electronics^ in part depends 
on the realization of devices such as molecular transis- 
tors, switches, or diodes. One strategy towards this 
goal involves the coupling of electronic and vibrational 
(phononic) degrees of freedom of molecules. Experiments 
with single-molecule devices have demonstrated effects of 
electron-phonon coupling in current-voltage characteris- 
tics (JVs)j2i2ii and a number of theoretical studies have 
investigated such features in JVs^^i^^ shot noise^**^ 
and the thermopower— as well as applications such as 
diodes^ and switches^ 

A question of principal importance for single-molecule 
devices are the consequences of nonequilibrium effects at 
finite bias. Strong nonequilibrium molecular vibrations 
can be beneficial in molecular devices, e.g., by enhanc- 
ing switching rates between molecular conformations. In 
other instances, they may hinder the operation of de- 
vices, in the extreme case by inducing dissociation of 
the molecule. Recent theoretical work shows that even 
within simple models, vibrational nonequilibrium has im- 
portant effects on IVs and shot noise^^*^ may induce a 
shuttling instability?i£*i& or lead to current flow charac- 
terized by a self-similar hierarchy of avalanches of large 
numbers of transferred electrons^ 

Recent numerical results by Mitra et al£ suggest 
that intriguingly, vibrational nonequilibrium becomes 
stronger as the electron-phonon coupling A decreases. 
Characterizing the vibrational nonequilibrium by the 
probability distribution of phonon excitations, these au- 
thors observe that the width of this distribution grows 
with decreasing coupling A. These numerical results are 
obtained within a minimal model describing transport 
through one molecular orbital, which is coupled to a sin- 
gle vibrational mode. 

In this paper, we first clarify the underlying mechanism 
for this nonequilibrium effect by developing an analytical 
theory. Our approach relies on a mapping to a Fokker- 
Planck equation, which becomes exact in the limit of 
weak electron-phonon coupling. This mapping predicts 



that the width at half maximum (WHM) of the phonon 
distribution diverges as A — > 0. Remarkably, the WHM is 
shown to scale as A _Q with bias-dependent, non-integer 
exponents a > 0. We confirm our analytical results by 
numerical simulations. 

"Real" single-molecule devices will typically involve 
features such as several vibrational modes, anharmonic 
vibrational potentials, and direct vibrational relaxation 
(e.g. due to radiation or interaction with the substrate), 
which are not fully captured by the minimal model. We 
therefore discuss how various such extensions of the min- 
imal model, which may be important for an accurate de- 
scription of experimental systems, modify our analytical 
findings. Specifically, we include anharmonic vibrations 
within a Morse-potential model which allows us to dis- 
cuss current-induced dissociation of the molecule. In this 
context, we show that the current-induced dissociation 
rate is governed by an interplay of the above-mentioned 
divergence of the width of the phonon distribution and 
a slowing down of the diffusion in phonon space as A 
decreases. 

The outline of the paper is as follows: In Sec. [n] we dis- 
cuss the nonequilibrium effects of weak electron-phonon 
coupling within the minimal model. The model is spec- 
ified in III Al and the resulting properties of the phonon 
distribution are derived in III Bl Nonequilibrium proper- 
ties of "real" molecules are discussed in Sec. lllll bv going 
beyond the minimal model. In particular, we address 
the effects of vibrational relaxation and the presence of 
several phonon modes as well as the situation of anhar- 
monic potentials in IIII Al IIII Bl and IIII CI respectively. 
Our conclusions are summarized in Sec. IIVI 



II. WEAK ELECTRON-PHONON COUPLING 
WITHIN THE MINIMAL MODEL 

A. Reduced model 

We investigate the nonequilibrium vibrational proper- 
ties of a molecule coupled to metallic source and drain 
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electrodes under finite bias within the following mini- 
mal modeliSiiiiiLiS Electronic transport is taken to result 
from sequential tunneling through one spin-degenerate 
molecular orbital with energy s, which is measured rel- 
ative to the zero-bias Fermi energies of the leads and 
which can be tuned by a gate voltage. As a minimal 
model, we consider a single vibrational mode with fre- 
quency lo. [Typical vibrational energies in molecules are 
of the order of 0.1 eV.] The system's Hamiltonian reads 

H = H moi + i/l ca ds 



-flmix, where 



U 



-Hmoi =en d + — n d {n d - 1) 

+ Xhiu{tf + b)n d + hwitfb + 1/2) 



(1) 



describes the molecular degrees of freedom, -f/i ea ds a free 
electron gas in the leads a — L, R (with creation opera- 



tors 4 PCT ), and 



,=L,R p, a 



{taC apa d 



h.c. 



(2) 



the tunneling between leads and molecule. 

We focus on the regime of strong Coulomb block- 
ade, appropriate when voltage and temperature are small 
compared to the charging energy U. The operator d a 
(dj.) annihilates (creates) an electron with spin projec- 
tion a on the molecule, n d = ^ a d\d a denotes the corre- 
sponding occupation- number operator. Vibrational exci- 
tations are annihilated (created) by b (b'). The electron- 
phonon coupling term can be eliminated by a canonical 
transformation, 2iA£ leading to a renormalization of the 
parameters e and U, and of the lead-molecule coupling 
t a — > t a cxp[— X(b> — b)]. From now on, we refer to the 
renormalized parameters as e and U. 

The coupling between molecule and leads is parameter- 
ized by the tunneling matrix elements and i^, and it is 
assumed to be weak in the sense that the energy broad- 
ening r of molecular levels is small, i.e. T <C fcsT 1 , Hu>, 
so that a perturbative treatment for _ff m i x in the frame- 
work of rate equations, as introduced in the context of 
Coulomb blockade phenomena^ is appropriate. We fo- 
cus on temperatures fcgTC Sw (corresponding to typical 
low-temperature experiments, see e.g. Ref. 3|). For sim- 
plicity, we assume a symmetric device with tjr, = £r = to 
and identical voltage drops of V/2 across each junction. 29 

Then, the occupation probability P™ for the molecular 
state \n,q) with n electrons and q phonons is determined 
by the rate equations 



dPl 
dt 



E 

n' ,q' 



\p n !w n ! 

J* q i qi- 



- p n w n ~ 

Q 1— 



(3) 



(The discussion of direct phonon relaxation is deferred 
until later in this paper.) The transition rates 
obtained by Fermi's golden rule are proportional to the 
square of the Franck-Condon (FC) matrix element 



dx cj)qi (x)<f>q (x — V2\e o sc ) , (4) 



i.e. the overlap of two harmonic oscillator wavefunc- 
tions <f) q {ps), shifted relative to each other by a dis- 
tance \/2A£osc, with electron-phonon coupling strength 
A, and vibrational oscillator length ^ osc = {K/Muj) 1 / 2 . 
Based on spectral data for diatomic molecules^ one can 
find electron-phonon coupling strengths ranging between 
A ~ 0.01 (BeO) and A ~ 5.4 (Kr 2 ). The FC matrix 
elements can be expressed as 

Mq^q, = {q\/Q\) 1/2 \ Q -^e- x2 l 2 L^ (A 2 ) sga(q-Qy- Q , 

(5) 

see e.g. Ref. [ll], where L^(x) denotes the generalized 
Laguerre polynomial. 



B. Phonon distributions for weak electron-phonon 
coupling 



For A« 1, Eq. J5J leads to 



\M, 



Q\ X 2A " 

q\ (AqlY 



(6) 



valid for q\ 2 , AqX 2 -C 1, where Q = ma,x{qi,q 2 }, q = 
minjqi, 52}, and Aq = Q — q. Therefore, the FC matrix 
elements and the transition rates W™^^ 1 decay rapidly 
with increasing Aq. Consequently, the vibrational state 
of the molecule is predominantly changed by processes 
for which q — > q' = q ± 1, and Aq = 1. Neglecting 
all other processes, Eq. describes a random walk in 
the space of phonon states q with g-dependent nearest- 
neighbor hopping rates. In this approximation, the ran- 
dom walker would eventually escape to infinity, as the 
rates for q ^ q + 1 are equal and grow with q. This 
implies that there is no steady-state phonon distribution 
P q = ^2 n Pq within this random-walk model. 

To derive the actual steady-state phonon distribution, 
it is therefore imperative to go beyond the random-walk 
model by including higher-order processes with Aq > 1. 
These may favor vibrational de-excitation processes since 
the applied voltage sets an upper limit to the increase 
(but not to the decrease!) in the vibrational excita- 
tion q by a tunneling event. For example, for e = 
the full voltage drop eV/2 per sequential-tunneling event 
can be converted into vibrational energy. Thus, Aq a = 
\eV/2fojj\ + 1 is the leading-order asymmetric process for 
which only de-excitation processes are permitted. (Here, 
[rj denotes the largest integer smaller or equal to r.) 

We can now derive the scaling of the steady-state 
phonon distribution P q with electron-phonon coupling 
A by balancing the diffusion process due to tunneling 
events with Aq = 1 [with diffusion constant ~ qA 2 , see 
Eq. ©] and the leading asymmetric drift process [with 
rate (qA 2 ) A?a , see Eq. JJjJ]. This leads to the balance 
equation 



qA 2 P„ 



(qA 



2\\q a p I 
r 1 



(7) 
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FIG. 1: (a) Scaling behavior of the phonon distribution width 
as a function of electron-phonon coupling for representative 
values of (eV/huj, e/Jko) equal to a: (3,0); 6: (5, 1); c: (5,0); 
d: (7, 1). (b) Width of the phonon distribution as a function 
of bias and gate voltage at A = 0.15 and ksT = 0.005/icj. The 
scaling of the width differs for the plateaus according to A: 
go ~ A -1 , B: qo ~ A~ 4,/3 , and C: go ~ A~ 



-3/2 



which implies a scaling law for the width qo of P q , namely 



a = 2(Aq a - 1)/Aq a . 



(8) 



This power-law scaling is nicely confirmed by numeri- 
cal results for P q as shown in Fig. Qfa). Remarkably 
for ksT <C fiLo the discrete dependence of the leading 
asymmetric process on bias and gate voltage implies fi- 
nite regions in the (V, e)-plane characterized by certain 
non-integer exponents a. This "phase diagram" is shown 
in Fig. ^b) where the wedge-shaped regions A, B, and 
C correspond to a = 1,4/3, and 3/2, respectively^ 

For the diamond-shaped regions along the line e = in 
Fig. [Jb), we can go beyond the derivation of this scal- 
ing behavior and obtain analytical results for the entire 
phonon distribution P q by a mapping to a Fokker-Planck 
equation. The derivation exploits the crucial observa- 
tion that for U — > oo and fc^T <C hu>, the transition 
s"^" w q -> q i factorize into a spin factor 



(1 + 8 n t t i)8\ n — „'i i and a phonon factor 

2 



'0 



= r„ - \M q ^ q ,\* [6{q+Aq a -l-q')+6(q-q'-Aq a )], 



(9) 

where t 1 = T/H = 2irp\to\ jh. In the stationary case, 
this implies the factorization P™ = P n P q , which allows 
us to derive the purely phononic rate equation 



= dPq/dt 



i' 



PqW q - 



(10) 



Since the phonon distribution becomes wide, we can take 
q to be continuous, expand P q > around q' = q up to sec- 
ond order, and keep only the leading-order contributions 
to diffusion and drift. In this way we obtain the Fokker- 
Planck equation 

with diffusion coefficient D(q) = 2qX 2 /r , drift coeffi- 
cient A(q) = [\ 2 ~c{q\ 2 ) A ^}/T , and c = 2A 9Q (A 9a !)- 2 . 
Remarkably, the stationary Fokker-Planck equation Hll|) 
can be solved analytically for any Aq a by the scaling 
ansatz P q = aX a f(X a q) with a normalization constant 
a. The universal function / is uniquely determined by 
Eq. (|11|) together with the boundary conditions /(0) = 1, 
/'(0) = 0, and we find 



/Or) = exp [~x A ^/b] 



(12) 



rates W, 



q~^q 



with b = i(Ag a !) 2 . Note that in particular, this analyti- 
cal result confirms the power-law scaling (JSJ) of the width 
of the phonon distribution. The power-law scaling (JSJ 
together with the analytical phonon distributions ljl2|l 
constitute the central results of this paper. The phonon 
distributions are nicely confirmed by numerical solutions 
of the full rate equations as shown in Fig. [5] 31 

In fully developed nonequilibrium, the width of the 
phonon distribution diverges with decreasing electron- 
phonon coupling A, and the resulting phonon distribu- 
tions are non- analytic in A. An important theoretical im- 
plication of this result is that in fully developed nonequi- 
librium, perturbation theory in the electron-phonon cou- 
pling parameter A is inadequate. Indeed, we find below 
that the radius of convergence of such an expansion in A 
would involve the direct vibrational relaxation rate. 



III. IMPLICATIONS FOR "REAL" MOLECULES 

Transport through "real" molecules will typically in- 
volve physics that goes beyond the minimal model. In 
particular, we will discuss (i) direct vibrational relax- 
ation, (ii) the presence of more than one vibrational 
mode, and (iii) anharmonic vibrational potentials. We 
show that while the exact scaling results for phonon dis- 
tributions are specific to the minimal model, nonequi- 
librium effects at weak electron-phonon coupling persist 
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FIG. 2: (Color online) Phonon distributions P q at bias volt- 
ages (a) eV = 3huu and (b) eV = 5hu>, plotted for several cou- 
pling strengths A at e = 0, ksT = 0.05ftw. The insets show 
that these distributions (approximately) collapse to universal 
curves given by the red (dark gray) smooth curves, which are 
the solutions of the the Fokker-Planck equation 111! specific 
to each voltage range nhuj < eV/2 < (n + l)hu. 



creases. As A decreases, we therefore expect that there 
exists a crossover coupling Ao: For A 3> Ao, the vibra- 
tional diffusion is limited by the drift in phonon space 
induced by the asymmetry between vibrational excita- 
tion and de-excitation, as discussed above. By contrast, 
for A -C Ao, the dominant limiting process is direct vi- 
brational relaxation, leading to a decrease of the width of 
the phonon distribution. This expectation is confirmed 
by numerical results as seen in Fig.|2a) which shows the 
width of the phonon distribution as a function of A for 
various relaxation rates. Fig.[3fb) shows the dependence 
of Ao on relaxation time r for Aq a = 1. Note that Ao 
grows only very slowly with increasing relaxation. While 
the Ao vs. t dependence is close to a power law with an 
exponent 1/4 - 1/3, no simple scaling can be expected. 
The reason is that the scaling suggested by the rate equa- 
tion amended by the relaxation term is incompati- 
ble with the scaling implied by the boundary condition 
A 2 ( Q+1 )/'(0) = -t /tSci at q = 0, where S = 1 + P°. 



B. Additional vibrational modes 

Typical molecules have many vibrational modes of dif- 
ferent vibrational frequencies. We expect that the sce- 
nario of the minimal model is most relevant to molecules 
whose lowest frequency mode happens to be weakly cou- 
pled. As this mode becomes highly excited, it would 
start to mix with other (higher-frequency) modes. In the 
simplest approximation, we can account for such mode 
mixing as a channel of direct vibrational relaxation so 
that the discussion of the previous subsection applies. 
Indeed, due to this mixing, vibrational energy can be 
distributed among different modes, which will generally 
tend to decrease phonon occupation numbers, similar to 
the vibrational relaxation discussed above. Under these 
conditions, such a weakly coupled vibrational mode may 
provide an efficient pathway to "pump" higher-frequency 
vibrations. 



and can be understood within the phonon diffusion model 
as long as the vibrational relaxation rate remains small 
compared to I / e. 



A. Direct vibrational relaxation 

Direct phonon relaxation can be included within 
the relaxation-time approximation by adding 
— \\Pq — Pq q J2 q > Pq>] to tne r.h.s. of the rate equations 
Q. Here, P® q denotes the equilibrium phonon distri- 
bution, which can be approximated by P° q = 5 q ,o for 
kgT <C huj. 

To understand the effect of direct vibrational relax- 
ation on the phonon distribution P q , it is important to 
note that the diffusion and drift processes in phonon 
space slow down as the electron-phonon coupling A de- 



C. Morse potential and dissociation 

So far, our considerations were based on the harmonic 
approximation for the phonon potential. We argue how- 
ever that wide phonon distributions are not an artefact 
of this approximation, but also appear for more realistic, 
anharmonic potentials. As an example, we investigate 
the effect of weak electron-phonon coupling for the Morse 
potential 



V{x) = D 



-20{x-x o ) 



2e 



-f3(x-x ) 



(13) 



where D > denotes the dissociation energy, p the in- 
verse range, and Xq the potential minimum. 

The Morse potential^ accurately describes the vi- 
brations of diatomic molecules and allows us to study 
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FIG. 3: (a) Width of the phonon distribution as a function 
of electron-phonon coupling strength for several vibrational 
relaxation times r and eV = 3hui, ksT — 0.05hcu, e = 0. 
For A > 0.25 all curves show the approximate qo ~ A -1 scal- 
ing. Below a relaxation-rate dependent crossover-point Ao 
the WHM is strongly suppressed due to direct relaxation, (b) 
Crossover-point Ao vs. relaxation time r. 



current-induced molecular dissociation. 23 This phe- 
nomenon has been explored experimentally in STM ex- 
periments with molecules on metal surfaces, and suffi- 
ciently high currents have been reported to lead to fast 
dissociation of the absorbed molecules S While this sce- 
nario with strong coupling between the molecule and the 
metal surface is distinct from the regime addressed in 
this paper, we remark that, in principle, the deposition 
of molecules on passivated surfaces^ can be exploited to 
study the weak-coupling regime as well. 

In analogy to the harmonic oscillator model, we as- 
sume that the potential energy curves for the neutral 
and singly-charged molecule have the same shape (i.e. D 
and (3 are fixed), but are shifted with respect to each 
other by Ax = \/2\£ osc . Building on e.g. Ref. |2^, we 
derive Franck-Condon matrix elements for the Morse po- 
tential, analogous to Eq. J3J, between bound states as 
well as between bound and continuum states. Details of 
this calculation are referred to Appendix . 

Specifically, we study the current-induced dissociation 
rate of the molecule as function of electron-phonon cou- 
pling and bias voltage by Monte-Carlo simulation. As- 
suming low temperatures and switching on the voltage 
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FIG. 4: Mean dissociation rate as a function of electron- 
phonon coupling strength at ksT = 0.005?kj for several bias 
voltages as obtained by Monte Carlo simulation with a Morse 
potential containing 10 bound states. 



at t = 0, the molecule starts in the phonon ground state 
and then evolves in time due to the tunneling dynam- 
ics. Given that transitions from the continuum back to 
bound states are negligible, we obtain an average disso- 
ciation rate Tdiss by recording the times tdiss,i required 
for reaching the continuum and averaging over samples. 
(We note that calculations of mean first-passage times 
for the highest-lying bound level give compatible disso- 
ciation times.) Typical results for dissociation rates for 
weak electron-phonon couplings (without relaxation) are 
depicted in Fig. 0] 

The maximum in the dissociation rate Tdiss vs. A can 
be understood as a direct consequence of a competition 
between the broadening of the phonon distribution and 
the slowing down of diffusion in phonon space. As A 
decreases from values of the order of unity, Tdiss first 
increases. This reflects the concurrent increase in the 
width of the phonon distribution. Beyond the maximum, 
Tdiss decreases due to the slowing down of diffusion in 
phonon space. Finally, the dissociation rate increases 
with voltage because of the increased width of the phonon 
distribution (see Fig. and the possibility of multiple- 
phonon excitations within one tunneling event. 



IV. CONCLUSIONS 

We have studied the current-induced vibrational 
nonequilibrium in single-molecule devices and found that 
remarkably, the width of the nonequilibrium phonon dis- 
tribution increases with decreasing electron-phonon cou- 
pling. We have identified regions in the bias voltage- 
gate voltage plane in which the width of the phonon 
distribution exhibits power-law divergences with decreas- 
ing A, with voltage-dependent non-integer exponents. In 
some representative cases, we are able to derive ana- 
lytical phonon distributions by a mapping to a Fokker- 



6 



Planck equation. These striking effects of current- 
induced noncquilibrium are found to have important im- 
plications in more realistic models which include direct 
vibrational relaxation and anharmonic potential surfaces. 
A very important conclusion from our work is that ap- 
proaches which are perturbative in the electron-phonon 
coupling A have to be assessed with extreme care in fully- 
developed nonequilibrium. Finally, we remark that re- 
cent experiments 26 show that the vibrational relaxation 
time can be as large as 10ns, in which case current- 
induced vibrational nonequilibrium becomes important 
for currents as small as lOpA. 
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and \J\ [which is the integer closest to and smaller 
than j] gives the number of bound states. This set of 
functions has three appealing properties*^ (i) It forms 
a discrete complete orthonormal basis enumerated by 
n = 0, 1, 2, . . .. (ii) The first \J\ + 1 functions form a ba- 
sis for the bound eigenstates of the Morse potential, all 
remaining functions span the space of continuum eigen- 
states. (hi) With respect to this basis, the Hamiltonian 
takes a particularly simple tridiagonal form. 

Denoting the bound and continuum eigenstates of the 
Morse potential by | ip q ) and | ipE ), respectively, we can 
now calculate the relevant FC matrix element for a tran- 
sition from a bound state | tp q ) into a continuum state 
\iPe) by 
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APPENDIX A: FC MATRIX ELEMENTS FOR 
THE MORSE POTENTIAL 

The calculation of dissociation rates within the Morse- 
potential model requires not only the determination of 
FC matrix elements between bound states [which is 
straightforward and can be found, e.g., in Ref. |27| . 
but also matrix elements between bound and continuum 
states. Although analytical expressions for the contin- 
uum eigenfunctions of the Morse potential are known^ 
their structure involves confluent hypergeometric func- 
tions with complex parameters, rendering a direct nu- 
merical evaluation of the FC matrix elements difficult. 

Instead, we make use of the c omp lete set of orthonor- 
mal functions introduced in Ref. 25, 



(3n\ 



T{2a + n) 



L 2 r r 1 (y)y a e-y/ 2 . (Al) 



Here, we denote y = (2j + l)e & x and a = j — [j\ . j is 
fixed by the Morse potential parameters to 



2j + 1 = y/8fiD/(3 2 h 2 , 



(A2) 



M, 



'Iq^E 



in 

= E 



(A3) 



oo 



a *Enfikm 



m=0n=[jj+l 



Here, the expansion coefficients a£„ and (3k m for contin- 
uum and bound eigenstates with respect to the {| <f> n )} 
basis are obtained through numerical diagonalization of 
the Hamiltonian. Finally, the FC matrix for the {| <p n )} 
basis are given by 



(-1)" {Aa) a T{n + n' + 2a) 

y / T(2a + n)r(2cr + n')n\n'\ 
(A4) 



(a - 1)"+"' 

(a+ l)n+n'+2a' 



F x (-nf, -n; -n - ri - 2a + 1 



(l-a)- 



where a = exp — -4= and 2F1 denotes the (Gaussian) 

V A+ 1 

hypergeometric function. In practice, one introduces a 
cutoff for the basis {| 4> n )}> leading to a discretization of 
continuum eigenstates. In order to ensure a sufficiently 
dense spacing of the spectrum close to the dissociation 
limit, we take into account ~ 10,000 basis functions. 
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